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ABSTRACT 

We study the interplay among cooling, heating, conduction, and magnetic fields in 
gravitationally-stratified plasmas using simplified, plane-parallel numerical simulations. 
Since the physical heating mechanism remains uncertain in massive halos such as groups 
or clusters, we adopt a simple, phenomenological prescription which enforces global 
thermal equilibrium and prevents a cooling-flow. The plasma remains susceptible to 
local thermal instability, however, and cooling drives an inward flow of material. For 
physically plausible heating mechanisms in clusters, the thermal stability of the plasma 
is independent of its convective stability. We find that the ratio of the cooling timescale 
to the dynamical timescale t coo \/t^ controls the non-linear evolution and saturation of 
the thermal instability: when t coo \/t& < 1, the plasma develops extended multi-phase 
structure, whereas when t coo \/ts > 1 it does not. (In a companion paper, we show that 
the criterion for thermal instability in a more realistic, spherical potential is somewhat 
less stringent, t coo \/t^ < 10.) When thermal conduction is anisotropic with respect to 
the magnetic field, the criterion for multi-phase gas is essentially independent of the 
thermal conductivity of the plasma. Our criterion for local thermal instability to produce 
multi-phase structure is an extension of the cold vs. hot accretion modes in galaxy 
formation that applies at all radii in hot halos, not just to the virial shock. We show 
that this criterion is consistent with data on multi-phase gas in the ACCEPT sample of 
clusters; in addition, when t coo \/t& > 1, the net cooling rate to low temperatures and 
the mass flux to small radii are suppressed enough relative to models without heating to 
be qualitatively consistent with star formation rates and x-ray line emission in groups 
and clusters. 

Key words: galaxies: evolution, galaxies: halos, galaxies: clusters: intracluster medium, 
ISM: kinematics and dynamics 



1 INTRODUCTION 

While the formation of dark matter halos can be understood 
via gravitational interactions alone, the combined effects of 
cooling and gravity are essential to galaxy formation (Rees & 
Ostriker 1977; Silk 1977; White & Rees 1978). This interplay 
remains poorly understood, however, because the dense plasma 
in many high-mass halos is predicted to cool and accrete far 
more rapidly than is observed (Peterson & Fabian 2006). As a 
result, theoretical models and numerical simulations routinely 
over-predict the amount of cooling and star formation in 
massive galaxies (e.g. Saro et al. 2006); this discrepancy is 
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an example of the well-known "cooling- flow problem." Some 
studies avoid the cooling-flow problem by focusing only on very 
hot halos (e. g. Sijacki & Springel 2006) or by "pre-heating" the 
gas to very high entropies (e. g. Oh & Benson 2003; McCarthy 
et al. 2004); this solution cannot work in general, however, 
because such hot systems are not representative of the cluster 
population (Cavagnolo et al. 2008). In particular, the central 
cooling time in many clusters is shorter than the Hubble 
time. Significant heating ("feedback") is required even at low 
redshift to suppress cooling in high-mass halos (Benson et al. 
2003), and thus to explain the observed cutoff in the galaxy 
luminosity function (Cole et al. 2001; Kochanek et al. 2001). 

Detailed x-ray observations of groups and clusters also 
highlight the need for significant heating of the intracluster 
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plasma. Though these objects contain large amounts of ra- 
diating plasma (Fabian 1994), their x-ray spectra indicate 
a paucity of material cooling below ~ 1/3 of the maximum 
temperature (Peterson & Fabian 2006). This demonstrates 
that most of the plasma radiates without actually cooling; 
i. e., an energy source heats the plasma at a rate similar to its 
cooling rate. 

Although heating dramatically suppresses cooling in 
groups and clusters, there is clear evidence for some cool 
gas in these systems. Studying this cold material can provide 
an important window into the heating mechanisms in groups 
and clusters and may help us understand how the balance 
between heating and cooling is maintained. The existence 
of a cold phase can be inferred from star formation (O'Dea 
et al. 2010), but it has also been directly imaged in a number 
of cases, revealing filamentary nebulae located tens of kilo- 
parsecs from the center of the potential (e. g. Fabian et al. 
2008; McDonald et al. 2010, 2011a). Despite more than five 
decades of study, the origin of these dramatic filaments has 
yet to be conclusively established: they have been interpreted 
as the remnant from an enormous, central explosion (Lynds 
& Sandage 1963; Lynds 1970), mass dropout from a cooling 
catastrophe (Fabian & Nulsen 1977; Cowie et al. 1980; Nulsen 
1986), debris from a high-speed merger of two gas-rich galaxies 
(Holtzman et al. 1992), or material dredged from the central 
galaxy by rising bubbles inflated by its AGN (Fabian et al. 
2003, 2008). Studies have shown, however, that the cold gas is 
highly correlated with short central cooling times in the hot 
intracluster plasma (e. g. Hu et al. 1985; Heckman et al. 1989; 
Cavagnolo et al. 2008; Rafferty et al. 2008, clearly illusrated 
in Voit et al. 2008), suggesting that its origin involves cooling 
of the intracluster medium (icm). 

In this paper, we investigate the possibility that the cold 
phase forms as a consequence of local thermal instability in a 
globally stable atmosphere. Though many authors (e. g. Fabian 
& Nulsen 1977; Nulsen 1986) have previously proposed that the 
filaments form via thermal instability, this idea has typically 
been analyzed in the context of a cooling-flow background. 
Subsequent analytic and numerical studies (e. g. Malagoli et al. 
1987; Balbus 1988; Balbus & Soker 1989; Hattori & Habe 
1990; Malagoli et al. 1990; Joung et al. 2011) showed, however, 
that the linear thermal instability is ineffective at amplifying 
perturbations in a cooling flow and concluded that it is unlikely 
unlikely to produce the cool filaments seen in many clusters. By 
contrast, the thermal instability is not suppressed in a globally 
stable atmosphere (Defouw 1970; Balbus 1986), which is now 
believed to be a better approximation to the thermal state 
of the icm. Quantitatively studying the thermal instability in 
this context has proven difficult because of the cooling-flow 
problem: studies that include cooling and gravity generally 
find that the plasma is globally thermally unstable, and that 
the entire cluster core collapses monolithically. 1 

We avoid the cooling-flow problem in this paper using a 



1 One-dimensional models of the ICM with simplified heating pre- 
scriptions can be stable or quasi-stable with episodes of heating and 
cooling (e.g. Guo & Oh 2008; Ciotti & Ostriker 2001); however, 
creating a realistic, stable model in two or three dimensions is 
significantly more challenging. Moreover, just as convection cannot 
be modeled in one dimension, the dynamics of cool, over-dense gas 
sinking through the hot atmosphere is absent in one-dimensional 
models. 



new strategy. Rather than attempting ab initio calculations of 
heating in clusters, we start from the observational fact that 
the ICM does not cool catastrophically. We therefore implement 
a phenomenological heating model that enforces approximate 
thermal equilibrium when averaged over large scales. We use 
this model to study the formation of multi-phase structure 
and we compare our results with archival data for groups and 
clusters. We find that the thermal stability of the plasma does 
not depend on its convective stability (see section 4.2). Instead, 
we find that the ratio of the thermal instability timescale t T i 
to the dynamical (or "free-fall" ) timescale tff governs the non- 
linear saturation of the local thermal instability: the plasma 
develops extended, multi-phase structure only where this ratio 
falls below a critical threshold (§5.1). This conclusion is not 
sensitive to significant perturbations about our idealized feed- 
back prescription (§5.4) and is unchanged even in the presence 
of very rapid thermal conduction (§7.3) (Though the threshold 
may depend somewhat on the geometry and initial conditions 
of the system; see section 8). 

This paper is the first in a series; here we present our 
model of local thermal instability and demonstrate its proper- 
ties and implications using linear theory and non-linear simu- 
lations. The aim of this paper is to develop an understanding 
of the essential physics of the problem and we therefore study 
stratified plasmas using idealized, plane-parallel calculations. 
In our companion paper (Sharma et al. 2012; hereafter Pa- 
per II), we present more realistic calculations of groups and 
clusters with spherical geometries and NFW halos. In both 
papers, we focus our analysis on the transition of material 
from the hot phase to the cold phase; we are not yet able to 
quantitatively predict any precise properties (such as sizes or 
luminosities) of the cold filaments produced via thermal insta- 
bility. We discuss in section 8 how our results can nonetheless 
be tested observationally. 

Because we put in by hand that hot halos are in approxi- 
mate global thermal equilibrium, our model provides no direct 
insight into how this balance is maintained. This is both a 
weakness and a strength of our current approach: though our 
setup is necessarily phenomenological, our results are not tied 
to any particular heating mechanism. Thus, we expect that 
our conclusions should apply to a wide range of systems, rang- 
ing in mass from galaxies to galaxy clusters. We return to 
this point in sections 5 and 8 and we study more physically 
motivated heating models in Paper II. Our present aim is not 
to identify a plausible solution to the cooling-flow problem, 
but rather to understand what implications a stabilizing heat 
source has for the local thermal stability and dynamics of the 

ICM. 

The structure of this paper is as follows. In section 2, we 
describe our model for the plasma, including our phenomeno- 
logical heating prescription. We describe our numerical method 
in section 3, linear theory results in section 4, and our primary 
numerical results in section 5. Section 6 provides a physical 
interpretation of the numerical results. For simplicity, we ini- 
tially ignore magnetic fields and thermal conduction in this 
paper; section 7 shows results including these effects. Finally, 
in section 8, we speculate on the astrophysical implications of 
our model and compare our results with observational data 
from the accept catalog (Cavagnolo et al. 2009). 
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2 PLASMA MODEL 

In this section, we describe our model for the cooling, heating 
and dynamics of the plasma in a dark matter halo. Due to the 
wealth of observations of the icm, we explicitly motivate our 
model for galaxy clusters, and some of the details we present 
in this section may not apply to galaxies. Nonetheless, our 
analysis is fairly general and we expect that some of our basic 
conclusions also hold massive galaxies (see Paper II for more 
details). 

We model the plasma as an ideal gas, sitting in the 
fixed gravitational potential of the halo and subject to both 
optically-thin radiative cooling and heating by a stabilizing 
feedback mechanism. In the interest of simplicity, we initially 
ignore both thermal conduction and the dynamical effect of 
the magnetic field; these effects are important in the ICM (see, 
e.g. McCourt et al. 2011 and references therein), but do not 
change our qualitative conclusions. We generalize our results 
to conducting, magnetized plasmas in section 7. 

The equations for the conservation of mass and momen- 
tum in the plasma, and for the evolution of its internal energy 



dp 

dt 



+ V.(pt;) = 0, 



dt 



(pv) + V • (pv ®v + PI) = pg, 



(la) 
(lb) 
(lc) 



where p is the mass density, v is the fluid velocity, denotes 
a tensor product, P is the pressure, I is the unit matrix, g is 
the gravitational field, T is the temperature, 



7-1 pmu 



In 



(2) 



is the entropy per unit mass, and d/dt = d/dt + v • V is the 
Lagrangian (or convective) time derivative. In equation 2, &b 
is Boltzmann's constant and pmu is the mean mass of the 
particles contributing to thermal pressure in the plasma. The 
functions % and C describe heating and cooling of the plasma, 
respectively; we explain our prescriptions for these processes 
in the following sections. 

2.1 Feedback 

The physical origin of heating in clusters remains uncertain, 
but it is simple to understand why our model requires the 
heating function 7i. Equation lc shows that the timescale for 
the icm to cool in a cluster is ~ nT \C — if H = 0, this 

timescale near the centers of many clusters can be orders of 
magnitude shorter than the Hubble time (this is the aforemen- 
tioned cooling-flow problem) . The continued existence of the 
icm in these clusters therefore strongly suggests that it is very 
nearly in thermal equilibrium, with % approximately equal 
to C when averaged over sufficient length- or time-scales. 2 
Nonetheless, the multi-phase structure seen in many clusters 

2 An alternative is a "cooling-flow" model, where the cooling gas 
flows inward and is replenished by continued accretion (see Fabian 
1994). This is not a viable alternative to heating, however, as these 
models over-predict the rate of gas cooling to low temperatures 
and the star formation rates in clusters by a factor of 10—1000; 
furthermore, the resulting density profiles are strongly disfavored 



(e.g. McDonald et al. 2010, 2011a) suggests that the thermal 
instability also operates. For the purposes of this paper, we call 
this behavior globally stable, but locally thermally unstable. 

The processes maintaining global thermal stability in 
clusters are not fully understood. The condition of global 
stability with local instability constrains the possible heating 
mechanisms, however, and suggests a phenomenological model 
for heating in the icm. This is a model in which 7i « C 
on average, but not 1~L — C identically. We adopt a specific 
implementation of this feedback model which simply fixes 
thermal equilibrium at all radii in our model halos. We set 



n = {c), 



(3) 



where (• • •) denotes a spatial average at a given radius. Thus, 
heating in our simplified model is a function only of r and t. 
By construction, this heating function ensures global thermal 
equilibrium at all radii in the plasma (precluding a cooling 
catastrophe), but permits the thermal instability to grow on 
smaller scales. It thus captures what we believe is the essential 
physics for the formation of multi-phase structure and meets 
our observationally-motivated requirements for the thermal 
stability of the icm. 

Equation 3 can be roughly motivated by positing a causal 
relationship between cooling on small scales and heating on 
large scales. Accretion onto a central AGN induced by cooling 
at larger radii is a promising mechanism for this "feedback" 
(Pizzolato & Soker 2005, 2010), and feedback from star for- 
mation could play a similar role in lower mass halos. Our 
specific heating implementation instantaneously balances cool- 
ing in every radial shell — in detail, this behaviour is non-local, 
acausal, and unphysical. Equation 3 is intended to mimic the 
end result of very effective feedback, but does not directly 
model the feedback process. Finding a physically-motivated 
heating mechanism that also leads to global stability is an 
important goal in the theory of the icm, but it is outside the 
scope of our present study. 

Our heating model is necessarily idealized, and it is impor- 
tant to separate tautological results (put in by hand) from the 
results which more generally reflect the global stability and 
local instability of the plasma. Though the subtleties of feed- 
back are likely to strongly affect the evolution of the plasma, 
we find that our qualitative conclusions are not sensitive to 
the precise form of our heating function. We demonstrate this 
in section 5.4 by applying spatial and temporal variations to 
equation 3. Moreover, the simulations in Paper II reach simi- 
lar conclusions using a very different setup. Thus we believe 
that the results derived using equation 3 capture some of the 
essential (and robust) dynamics of local thermal instability in 
globally stable systems. It is, however, difficult to prove this 
conclusively given current uncertainties in the heating of the 

ICM. 

In our heating model, spatial variations between heating 
and cooling drive thermal instability; a more realistic model 
would likely introduce temporal, in addition to spatial, vari- 
ations. We show in §5.4 that our conclusions do not change 
unless these temporal fluctuations around the thermal equilib- 
rium are very large (~ 300%). We choose to begin our study 
using equation 3 because it is analytically tractable and lends 
itself to a thorough investigation. 

by x-ray observations (McNamara & Nulsen 2007). We therefore 
do not consider cooling-flow models in this work. 
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Equation 3 is appropriate for a heating process which dis- 
tributes energy per unit volume, such as the dissipation of MHD 
waves. Other processes like photoelectric heating distribute 
energy per unit mass. Since it is not yet known how feedback 
energy is thermalized in the icm, we generalize equation 3 to 
other processes: 



H = n 



<n«>" 



(4) 



Here, a = corresponds to volumetric heating and a = 1 
corresponds to mass- weighted heating. We show in section 4.1 
that the thermal instability takes ~ 3 times longer to develop 
in plasmas with a = 1 than in plasmas with a = 0; after scaling 
the timescales by this factor, however, we find very similar 
evolution for plasmas with volumetric and mass-weighted 
heating (see Fig. 2, below). 



2.2 Cooling 

In the idealized spirit of this paper, we adopt a simple cooling 
function C dominated by thermal Bremsstrahlung 



C B = n A(T) = A n T 



1/2 



(5) 



where n = p//nmn is the number density of particles in the 
plasma and we have introduced the standard notation A(T) for 
consistency with other work. Our conclusions are not sensitive 
to the shape of the cooling function as long as the plasma 
remains locally thermally unstable (§4); this is the case in the 
icm for temperatures above ~ 10 4 K. 

In an unstratified plasma, thermally unstable clumps of 
cool gas collapse to the Field length in the cold phase (the 
length-scale below which thermal conduction suppresses local 
thermal instability; Field 1965). Resolving the realistic Field 
length in the cold phase of the icm is numerically impractical, 
so we introduce a temperature floor at which we truncate the 
cooling function (see Sharma et al. 2010, §2.2 for a discussion 
of this approximation; also see §5.3 of this paper). We use the 
modified cooling function 



C — Cb Oh(T — Tfloor), 



(6) 



where Bh is the Heaviside function, and Tfl oor effectively 
becomes the temperature of the cold phase. 

The microphysical processes heating and cooling the cold 
phase in the icm are likely to be very complicated (see Ferland 
et al. 2009) and we do not consider them here. Our use of 
a temperature floor amounts to the reasonable assumption 
that, once a thermally-unstable fluid element cools below 
Tfloor, it is unlikely to enter back into the hot phase. This 
simplification, along with our omission of line-cooling from 
equation 5, prevents us from studying the evolution of the cold 
material in detail. This is not a major limitation, however, 
because we are primarily interested in the transition of material 
from the hot phase to the cold phase. Following the internal 
structure of the cold clumps would be crucial for calculating 
the emission from filaments or for studying the intermittency 
in the accreted mass flux, but these applications are beyond 
the scope of our present study. 

The simplified cooling function used here (eq. 6) prevents 
the gas from cooling below Tfloor and therefore artificially low- 
ers the gas density in thermally unstable clumps or filaments. 
We have confirmed, however, that the quantitative results in 



this paper are insensitive to the numeric value of Ta 

oor; pro- 
vided it is much lower than the initial (or virial) temperature 
of the plasma. In Paper II we use a realistic cooling function 
that includes both Bremsstrahlung and line emission, and 
which does not implement a temperature floor. The results 
from this more realistic model agree with our conclusions here. 



3 NUMERICAL MODEL 

We solve equations la-lc using the conservative MHD code 
Athena, modified to implement equations 4 and 6 via a 
semi-implicit, operator-split method (Sharma et al. 2010). 
Specifically, we evolve the thermal energy per unit volume 
E — n k&T/ (7 — 1) using 



SE M = (n (n) -£ (n) ) St 
< E (n) + SE (n) 



SE {n) > 



?; (n) /(i + |^ (n) /^ (n) l) ^ (n) 



< 



(7a 



(7b) 



where S indicates a finite approximation to a differential, and 
denotes the function / during the n th time-step of the 
simulation. This method explicitly prevents the temperature 
from becoming negative, even in the extreme case that the 
cooling time becomes shorter than the simulation time-step 
(although equation 7 is no longer accurate in this limit). Equa- 
tion 7 is asymmetric and is only accurate to first order in 
St /t cooi. In order to test the sensitivity of our simulations 
to these shortcomings, we have also run simulations using a 
fully explicit, sub-cycled method. The two methods yield very 
similar results. We use equation 7 because it is faster than an 
explicit method and because it does not alter our results. 

We perform most of our calculations on 2d Cartesian 
grids of resolution (300) 2 or 3d Cartesian grids of resolution 
(128) 3 . We show a resolution study in section 5.3. In the 
remaining sections, as in our simulations, we work in units 
with /cb = /im p = 1. 

We perform our calculations in the plane-parallel approxi- 
mation, with g — —g(z) z. We therefore use the words 'height' 
and 'radius' interchangeably in the following sections. We 
make our setup symmetric about the z = plane, with 



9 = 9o 



z/a 



[1 + {z/ay\ 



211/2' 



(8) 



Thus, g is nearly constant outside \z\ = a, with a smooth 
transition through zero at the center. This setup enables us to 
place the computational boundaries far from the center, where 
most of the cooling and feedback take place (see Fig. 1). To 
further diminish the influence of the boundaries, we end our 
simulations before one cooling time transpires at the boundary. 
We use reflecting boundary conditions in the direction parallel 
to gravity and periodic boundary conditions in the orthogonal 
directions. 

We set the softening radius a = 0.1 H, where H is the 
plasma scale-height (defined below). We turn off cooling and 
heating within \z\ < a because the physics at small radii is 
particularly uncertain and our feedback prescription (equa- 
tion 4) may not be a good approximation to what happens 
there. We allow cold material to accumulate in the center 
\z\ < a, but we otherwise ignore this region in our analysis. 
We have also performed simulations in which we do not turn 
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off cooling in the center and have confirmed that it does not 
change our conclusions at larger radii z > H. 

We initialize the ICM in hydrostatic equilibrium, with a 
constant temperature To and with the density profile 



p(z) = po exp [-J([l + (z/a) 2 ] 1/2 - l)] 



(9) 



where the scale-height H = To /go. For computational conve- 
nience, we set po = To = go = 1 and we take Tfl 00 r — 1/20. 
This roughly corresponds to a virialized halo, in which the ther- 
mal and gravitational energy in the plasma are approximately 
equal. The atmosphere defined by equation 9 is buoyantly 
stable, with ds/dz > 0. To test the sensitivity of our results to 
stratification, we also use the buoyantly neutral atmosphere 
defined by 



T(z) 



7 — 1 a 

7 

1/(7-1) 



|([l + (,/a) 2 ] 1/2 -l) 



(10a) 
(10b) 



We refer to the conditions defined by equations 9 and 10 
as isothermal and isentropic, respectively. Note that our use 
of the entropy gradient to determine convective stability is 
only appropriate because we have neglected conduction. When 
thermal conduction is efficient, the temperature gradient and 
the magnetic field orientation determine the convective stabil- 
ity of the plasma (Balbus 2000; Quataert 2008). We describe 
this in more detail in section 7. 

We seed thermal instability in our model atmospheres 
by applying an isobaric perturbation with a fiat spectrum 
ranging from k = 2n / L to k = 40tv/L, where L is the size 
of the simulation domain. The cutoff at high k makes the 
perturbation independent of resolution and permits a detailed 
convergence study. Unless otherwise noted, the modes of this 
perturbation have Gaussian-random amplitudes with an RMS 
value of 10 -2 . 

We define the free-fall time and the cooling time as follows: 



tB = ( — 

,30 



1/2 



2 nAn ' 



(11a) 
(lib) 



Since these timescales are functions of height in our simula- 
tions, we quote them in the plane z = H to give single values. 
When our analysis depends on the ratio t-n/ts, we restrict it 
to this plane. We perform simulations with different initial 
values of the ratio t coo \/tff by changing the parameter Ao; 
this permits direct and unambiguous comparison among our 
simulations because each is initialized identically. In reality, of 
course, Ao is set by fundamental physics, and different values 
of tcooi/tff correspond to clusters with different icm entropies 
or densities. 



4 LINEAR THEORY RESULTS 

Equations 1-6 completely specify our model. In the rest of 
this paper, we study the properties of this model and apply 
it to astrophysical systems. In this section, we describe the 
linear stability of our model and derive the timescale for the 
formation of multi-phase structure in the plasma. We discuss 



the well-known linear results in some detail because they 
inform our interpretation of the non-linear behavior described 
later. In addition, the interpretation of these linear results 
has generated some confusion in the literature, leading to 
conflicting claims about the thermal stability of gas in hot 
halos. 



4.1 Linear Stability 

We define the net cooling rate = C — Ti and assume that the 
plasma is initially in thermal equilibrium with = every- 
where. The derivative (dS/dT)p describes how the net cooling 
responds to a linear, Eulerian perturbation. If this derivative 
is negative, a decrease in temperature at a fixed location in 
the plasma leads to an increase in the net cooling rate; thus, 
the temperature decreases further and the perturbed fluid 
element runs away to low temperatures. Similarly, an increase 
in temperature causes the fluid element to run away to high 
temperatures. The plasma is therefore unstable to small tem- 
perature fluctuations when (<90/<9T)p < 0. A similar line of 
reasoning demonstrates that the plasma is thermally stable if 
(dS/dT) P > 0. Following Field (1965), we derive this result 
by linearizing and perturbing equations la— lc. This analysis 
yields the linear growth rate of the perturbations and will 
assist our interpretation of the non-linear results presented 
later. 

We Fourier transform equations la— lc and perform a 
standard WKB analysis. We seek solutions with growth times 
much longer than the sound-crossing time and therefore make 
the Boussinesq approximation, which filters out sound waves 
(Balbus 2000, 2001). Under these approximations, the dynam- 
ical equations become 



k • 6v = 



-iu)5s + Sv. 



iujk 5v = 
ds 



dz 



Sn 
n 

SB 
' nT 



[k 2 g-k(k-g)} 



(12a) 
(12b) 

(12c) 



In deriving equation 12b, we have crossed the momentum 
equation with k twice and used equation 12a to eliminate the 
compressive component of the velocity. This is consistent with 
the Boussinesq approximation and simplifies the algebra later 
on. Additionally, in the Boussinesq limit, 

Sa = -^u*i (13) 



7 — 1 n 



and 



£0 



dO\ 6n 
dTj p n' 



(14) 



In deriving equation 14, we have used the thermodynamic 
identity 



d In X 
<91nT 



d In X 
<91nT 



d In X 

dhin 



(15) 



for any state function X(n,T). Note that, although the net 
cooling rate varies explicitly with position, this dependence 
does not enter into equation 14 because £0 represents an Eule- 
rian perturbation at a fixed point in space. Although heating 
in our model has explicit radial and temporal dependencies, it 
experiences no first-order change under an Eulerian perturba- 
tion. We therefore ignore changes to the heating in this linear 
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analysis. This should not give the impression that heating is 
immaterial to the linear results; on the contrary, these results 
presume an initial equilibrium state with a stabilizing heat 
source. The growth of the thermal instability is very different 
in the absence of such heating (Balbus 1988; Balbus & Soker 
1989). 

Combining equations 12-14, we find that the linear dis- 
persion relation for the plasma is 

co 2 - i I ) c^cooi - N 2 (l - kl) = (16) 



where 



7 - 1 £> ( . x-i 

^cool — — — {'Jtcool) 



7 nT 



is the cooling rate, 



N 



— 1 ds 

9 

7 oz 



(17) 



(18) 



is the frequency for internal gravity waves, and k = k/k is 
the direction of the wave vector of the perturbation. As noted 
previously, we have neglected conduction and thus equations 16 
and 18 only apply on relatively large length scales, > the Field 
length (see section 7). Equation 16 implies that perturbations 
grow exponentially in amplitude ~ e PTlt , with 



7 " 



7 

7-1 

7 

3 _ 
2 



11 

n 



_ din A 
~ d\nT 



nT 



OL CJcool- 



(19a) 
(19b) 
(19c) 



The three forms of equation 19 are equivalent and are useful 
in different contexts. In equation 19c, we have specialized to 
Bremsstrahlung cooling. In this case, plasmas with a < 3/2 
are locally thermally unstable (with p T1 > 0), even though our 
model is (by construction) globally stable against a cooling 
catastrophe. 



4.2 Local Stability, Global Stability, and 
Convection 

Following Field (1965) and Defouw (1970), we showed in the 
previous section that the icm is likely to be locally thermally 
unstable, and we propose that thermal instability may produce 
at least some of the multi-phase structure in galaxy clusters. 
At first, our analysis may appear inconsistent with other 
claims (such as can be found in, e. g. Balbus & Soker 1989 
and Binney et al. 2009) about the importance of local thermal 
instability in galaxy and cluster halos. We review this apparent 
contradiction here and show that there is no inconsistency. 

Balbus (1988) and Balbus & Soker (1989) extensively 
studied thermal instability using Lagrangian techniques and 
discovered that it is significantly stabilized in a cooling-flow. 
In a globally stable atmosphere, however, perturbations do 
grow exponentially (Defouw 1970; Balbus 1986). Since it is 
now thought that clusters are globally thermally stable and 
that the icm persists for many cooling times, we expect the 
thermal instability to undergo many e-foldings and to become 
highly non-linear in clusters (though this does not always 
imply a large amplitude; see §6). Thus, assumptions about 



the global stability of the icm also dictate conclusions about 
its local thermal stability and one must be careful to choose 
an appropriate background model. 

Even though we expect perturbations to grow exponen- 
tially in clusters, they do not necessarily grow monotonically: 
equation 16 shows that a thermally unstable perturbation 
oscillates as it grows if the cooling time is longer than the 
buoyancy time. This overstability represents a driven gravity 
wave (see Defouw 1970). Since the thermal instability in this 
case is not purely condensational, its identification with multi- 
phase gas becomes somewhat unclear (Malagoli et al. 1987; 
Binney et al. 2009). However, the growth rate of the thermal 
instability is essentially unaffected by buoyancy (equation 16), 
and thus perturbations are also likely to become highly non- 
linear in this limit. Earlier studies of the thermal overstability 
in stratified plasmas have either focused entirely on the linear 
evolution of perturbations (Defouw 1970; Malagoli et al. 1987; 
Binney et al. 2009) or have studied them in the context of a 
cooling-flow (Hattori & Habe 1990; Malagoli et al. 1990; Joung 
et al. 2011), in which the thermal instability is suppressed 
(Balbus & Soker 1989). 

For the reasons listed above, we argue that earlier studies 
cannot directly predict the astrophysical implications of ther- 
mal instability in cluster halos. The astrophysical implications 
of the thermal instability depend on how the linear growth 
saturates in a globally stable environment. This motivates 
our present study. A series of previous investigations are very 
similar to ours (Nulsen 1986; Pizzolato & Soker 2005; Soker 
2006; Pizzolato & Soker 2010), but focus on the survival of 
preexisting cold filaments rather than their formation via ther- 
mal instability. Our investigation compliments these studies 
and produces the initial conditions they require. 

The saturation of the thermal instability involves the 
sinking of cool over-densities; in this respect, it bears some 
similarity to convection. This connection between thermal and 
convective stability was first recognized by Defouw (1970) and 
was significantly sharpened by Balbus & Soker (1989). Specif- 
ically, Balbus & Soker (1989) showed that thermal instability 
necessarily implies convective instability if the heating and 
cooling are state functions of the plasma. Heating in galaxy 
groups and clusters is very unlikely to be a state function of 
the ICM plasma, however. As a concrete example of spatially 
dependent heating, consider heating by turbulence (induced 
by, e. g. buoyant bubbles created by star formation or an agn). 
The heating rate in this case is set by the rate at which tur- 
bulent energy is transferred to small scales, and thus by the 
turbulence properties as a function of position. In this case, 
i. e. when the heating depends explicitly on position, there is 
no one-to-one relationship between convective and thermal 
stability (as noted by Balbus & Soker 1989). Fundamentally, 
buoyancy determines convective stability, while heating and 
cooling determine thermal stability; these processes are not re- 
lated in a globally stable atmosphere, and the thermal stability 
of an atmosphere is independent of its convective stability. 



5 SIMULATION RESULTS 

We extend our analysis into the non-linear regime using the 
numerical setup described in section 3. We have run a large 
suite of 2d and 3d simulations, summarized in Table 1. We fo- 
cus our analysis on the presence of multi-phase structure (§5.1) 
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Figure 1. Snapshots of the density (top) and fractional density inhomogeneity Sp/p = (p — (p))/(p) (bottom) at the time t = 10t T j(z = H) 
in our simulations. Note that our simulations are symmetric about the plane z = 0; this enables us to put the boundaries far from the center, 
where most of the cooling and feedback take place. Gravity points down in the top half of the domain and up in the bottom half (for the 
remainder of the paper, we primarily show images of the top of the domain). We applied the heating function in equation 4 to ensure global 
thermal stability, distributing the energy per unit volume (a = 0). From left to right, these simulations have initial values of t T i/% = 10, 3, 1 
and 0.1. These simulations demonstrate that cooling and heating drive internal gravity waves when tn/ts > 1. The amplitude of these 
waves increases with the cooling rate and approaches the size of the simulation domain when t TI ~ %. When £ti/% < 1, the thermally 
unstable gas collapses into dense clumps, which then rain down into the center of the potential. For clarity, we have restricted the color bar 
on plots with £ T i/tff = 1/10. In this simulation, (5p/p) max ~ 20 is set by the (arbitrary) temperature floor we impose. While the simulations 
have been run for 10 thermal instability times at z = if, gas near the boundaries has not yet had time to cool. The initial perturbations are 
still visible near the boundaries in the lower-rightmost plot. This figure also clearly shows the accumulation of cool material in the center 
of our simulation domain. We describe this process in more detail in §6. Animated versions of the figures in this paper can be found at: 
http:/ /astro. berkeley.edu/~mkmcc/research/thermal_instability/movies. html. 



Table 1. Parameters for simulations without conduction (§5). 



Initial 
Condition 


L/H 


a 




Isothermal 


3 





1.57 A" 1 




1 


4.70 A" 1 


Isentropic 


2 





1.25 A" 1 




1 


3.74A" 1 



We performed all simulations on square Cartesian grids of resolution 
(300) 2 or (128) 3 and physical size 2L (the scale-height H is defined 
in §3). We also performed simulations at other resolutions as part 
of a convergence study (§5.3). The cooling constant Ao is a free 
parameter in our model, which we choose to obtain the desired 
tn/tff. Each combination of the listed parameters was simulated 
with initial values of log 10 (£ T i/%) at z = H spanning between — 1 
and 1 with increments of 1/4. The top row represents our fiducial 
setup; we also performed 3d simulations using this setup with 
log 10 (_Ti/tff) = -1, -0.75, -0.5, 0, and 1. 

and on the accreted mass flux (§5.2), both of which can be 
compared with observations of groups and clusters. 

Equation 19c shows that the growth rate of the thermal 
instability is a factor of 3 smaller in plasmas with heating 
per unit mass than in plasmas with heating per unit volume. 
More generally, the timescale depends on the uncertain pa- 
rameter a and cannot be directly applied to (or inferred from) 



observations. Nonetheless, it is convenient to use t TI = 
to normalize time when considering the physics of the ther- 
mal instability with different values of a. We also use the 
cooling time t coo \ = E/ C — (7cj coo i) -1 when we compare our 
results with observations. These two timescales differ only by 
an uncertain factor of order unity. 

5.1 Multi-phase Structure 

We performed simulations with the ratio of time-scales £ T i/£ff 
ranging from 0.1 to 10 (measured at z — H) and ran each for 
ten growth times, until t = 10 t TI . Fig. 1 shows representative 
snapshots of the density at the end of our fiducial simulations 
with volumetric heating and isothermal initial conditions. Our 
simulations show that plasmas with cooling times shorter 
than the dynamical time (t T i <tff) develop spatially extended 
multi-phase structure, whereas plasmas with cooling times 
longer than the dynamical time (t TI > tff ) do not. Thus, the 
ratio tn/tft controls the non-linear saturation of the thermal 
instability in stratified plasmas; this observationally-testable 
prediction is the primary result of our study. This conclusion 
does not depend strongly on either the initial stratification of 
the plasma or on our choice of heating per unit volume. To 
demonstrate this, Fig. 2 shows variations of our fiducial simula- 
tions with isentropic initial conditions and with mass- weighted 
heating. In all four cases, the saturated state transitions from 
single-phase to multi-phase when the ratio of time-scales t T i / tff 
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Figure 2. Snapshots of the density at the time t = 10f T i, for different values of the time-scale ratio t T i/%- The top two rows show our 
simulations with isothermal initial conditions and volume and mass-weighted heating, while the bottom two rows show isentropic initial 
conditions with volume and mass-weighted heating. These results show that the non-linear behavior of the thermal instability is relatively 
independent of the initial stratification and the details of the heating. Note that the color scale was chosen to show the features in the gas 
and varies from plot to plot; the ranges for a given value of tri/% are similar to those shown in Fig. 1. 



becomes less than one. Below, we describe the plasma prop- 
erties in these two limits and the physics of the transition 
between them. 

The evolution of plasmas with short cooling times t Ti <C tff 
is straightforward: in this limit, the thermal instability devel- 
ops and saturates before the plasma can buoyantly respond. 
The initial perturbations therefore collapse into dense clumps 
essentially in-situ, and the ICM develops a highly inhomoge- 
neous, multi-phase structure wherever t Ti (z) < t. The clumps 
of cold gas then rain down onto the central galaxy on the 
(much longer) free-fall time, while bubbles of heated gas rise 
outwards. The rightmost panels of Figs. 1 and 2 illustrate 
this behavior. The result is a hotter atmosphere (in which 
tri/tff > 1), filled with clumps of cold gas. We show in Paper II 



that this end state resembles the observed properties of some 
cool-core groups and clusters. 

The saturation of the thermal instability is fundamen- 
tally different when the cooling time is long compared to the 
dynamical time. In this limit, gravity and buoyancy influence 
the linear evolution of the perturbations (though the growth 
rate changes only by a factor of two). Nonlinearly, however, 
buoyancy provides a critical saturation channel for the thermal 
instability that prevents the formation of multi-phase gas. This 
conclusion is qualitatively similar to that reached by Balbus 
& Soker (1989); however, the physics is very different in our 
case because the background atmosphere remains statistically 
in thermal equilibrium for many cooling timescales. As initial 
perturbations cool and grow, they sink in the gravitational 
potential and mix with gas at lower radii. The cooling thus 
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tn/tff *Ti/*ff 

Figure 3. (Left:) Mass fraction of cold material (with T < Ti n i t i a i/3) as a function of the timescale ratio tn/ts- The mass in cold material 
drops off sharply when t TI ~ %, and there is no extended multi-phase structure in the weak-cooling limit. All quantities in these plots 
represent averages from z = 0.9-1.1 if and from t = 9-10 t T i- (Right:) Fractional density inhomogeneity Sp/p as a function of the timescale 
ratio tri/%- Blue (red) lines indicate isothermal (isentropic) initial conditions, solid (dashed) lines indicate volumetric (mass- weighted) 
heating. The blue stars represent 3d simulations using our fiducial setup (isothermal initial condition with volumetric heating); the remaining 
simulations are 2d. 



drives a slow, inward flow of material; the associated mass flux 
is, however, significantly smaller than is predicted by models 
without heating. We return to this point in the following sec- 
tions. Rather than creating strong density inhomogeneities, 
cooling in this limit excites internal gravity waves with an 
amplitude that depends on the timescale ratio tn/ts- These 
waves represent the overstability highlighted by Balbus & 
Soker (1989) and by Binney et al. (2009); we discuss their 
saturation below. 

Our results depend crucially on the existence of a globally- 
stabilizing heating mechanism; if heating were not present, the 
atmospheres shown in Figs. 1 and 2 would collapse monolithi- 
cally. This globally unstable case has been studied extensively 
by Balbus & Soker (1989). Consistent with their analysis, 
we find that atmospheres with small initial density inhomo- 
geneities do not form multi-phase gas, regardless of tn/ts (see 
Paper II for a more detailed discussion). 

Figs. 1 and 2 show that, assuming the existence of a 
globally stabilizing heating mechanism, plasmas with short 
cooling times t Ti <C ts develop spatially extended multi-phase 
structure, while plasmas with long cooling times do not. The 
left panel of Fig. 3 demonstrates this result more quantitatively. 
Here, we plot the mass fraction of cold gas (with T < 1/3 To) 
at late times in the plane z — H as a function of t TI / tff . (Recall 
that, in our units, To ~ T v i r iai.) This figure shows that the 
fraction of cold gas drops precipitously around tn/ts ~ 1 and 
that there is essentially no multi-phase gas at large radii in 
simulations with tn/ts > 1. 

The right panel of Fig. 3 quantifies the dependence of the 
saturated density fluctuations on the time-scale ratio tn/ts 
and hints at the physics of the transition between the two 
limits. Here, we plot the root-mean-square (rms) average of 



the density perturbations 
Sp 

p 



(p) 



(p) 



(20) 



as a function of tn/ts in the plane z = H as in equation 3, 
(• • •) indicates a spatial average at a given radius. In the short 
cooling time limit, the plasma develops multi-phase structure 
with large density perturbations Sp/p > 1. By contrast, in the 
long cooling time limit, the density perturbations saturate at 
much lower values Sp/p <C 1. For plasmas with stable back- 
ground stratification (e.g. our isothermal initial conditions), 
Sp/p in this limit represents the amplitude of the gravity waves 
driven by cooling. Note that, while the mass fraction of cold 
gas drops off sharply around tn/ts ~ 1, Fig. 3 shows that the 
mean density fluctuation is a smooth function of this parame- 
ter, even in the weak cooling limit. We work to understand 
this quantitatively in section 6. 

We emphasize that the difference in the evolution of 
plasmas with long and short cooling times does not simply 
result from the thermal instability taking longer to develop in 
simulations with weak cooling. We have run each simulation 
for a fixed number of growth times t TI and, if gravity were 
not present, the results of our simulations with rapid and 
slow cooling would be nearly identical (we have confirmed 
this numerically). In fact, even in our simulations with long 
cooling times, the density contrast Sp/p becomes large near 
the center of the potential, where gravity is weak (eq. 8) and 
the thermal instability has time to develop. The development 
of multi-phase structure depends on both gravity and cooling 
and therefore on the parameter tn/ts, rather than simply on 
the cooling time alone. 
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Figure 4. Mass flux (averaged from z = 0.9-1.1 H) as a function of 
time in 3d simulations, normalized to the cooling-flow flux for that 
atmosphere. The mass flux is severely suppressed when t TI > % (see 
eq. 27). The mass flux is not suppressed as strongly when t TI < %, 
but it is highly variable. Thick lines indicate a positive mass flux 
(i. e., an outflow), while thin lines indicate a negative mass flux (i. e., 
an inflow). Note that gravity waves dominate the instantaneous 
mass flux when tn/ts = 10; the time-averaged accretion rate is 
much smaller than suggested by this plot. 



5.2 Accreted Mass Flux 

Fig. 4 shows the instantaneous, mean mass flux through the 
plane z — H as a function of time in three of our 3d, fiducial 
simulations. The mass fluxes are normalized to the values pre- 
dicted by cooling- flow models without heating, Mcf = pH/tn. 
It is clear that the mass flux is strongly suppressed relative 
to the cooling- flow solution whenever tn/ts > 1. This is not 
a trivial consequence of our feedback heating mechanism, be- 
cause for £ TI < tft , M approaches Mcf- Rather, the suppression 
of M for tn > tfi is also due to the non-linear saturation of 
the thermal instability (described below). The mass fluxes we 
find for t TI > ts are < 1% of the cooling- flow estimates and 
are therefore reasonably consistent with observational limits 
for cooling in the icm (Peterson & Fabian 2006). In Paper II 
we show that this suppression is even stronger in spherical 
potentials and we explore its dependence on the details of our 
heating model. 

In the rapid cooling limit, we find that gas heated at small 
radii, where the cooling time is shorter, rises up through the 
plane z — H and initially drives an outflow. As the thermal 
instability progresses, however, this outflow reverses and a 
strong accretion flow develops (although the accreted mate- 
rial is all in the cold phase, rather than the hot phase; see 
the left panel of Fig. 3). The accretion rate approaches the 
cooling-flow value and eventually depletes the atmosphere 
of its gas. Thus, even our idealized feedback model (eq. 4) 
cannot suppress a cooling catastrophe when tn/ts <C 1. We 
discuss the implications of this result in section 8 and, more 
thoroughly, in Paper II. 

Fig. 4 is only meant to be suggestive, as several sub- 



tleties in our analysis complicate a precise interpretation of 
the accreted mass flux. For instance, we use only the initial 
value Mcf, though this quantity changes dramatically over 
the course of some of our simulations. Additionally, in simula- 
tions with t T j <C tff, a more appropriate normalization for the 
mass flux might be pH/ts, since the gas is not likely to flow 
in faster than its free-fall rate. Paper II presents a much more 
realistic and thorough analysis of mass accretion rates. 



5.3 Resolution Study 

We test the numerical convergence of our results with 2d calcu- 
lations on grids of resolution (100)? (200)? (300)? and (400) 2 
for the full range of tn/ts- Additionally, we performed 3d 
calculations of resolution (128) 3 and (256) 3 for atmospheres 
with tn/ts = 0.1, 1, and 10. Fig. 5 shows the distribution of 
density perturbations in our 2d simulations; this quantity has 
no apparent trend with resolution. Similarly, our 3d simula- 
tions are nearly identical at resolutions of (128) 3 and (256) 3 . 
Though Fig. 5 only demonstrates convergence of an integrated 
quantity, our simulations also "look" very similar at different 
resolution: for example, in the rapid-cooling limit, the clumps 
of cold gas have similar shapes and sizes, and they appear in 
the same locations. 

We find rapid convergence in our simulations, even with- 
out including thermal conduction. By contrast, Sharma et al. 
(2010) found that convergence requires resolving the Field 
length (§7.2) in the cold phase of the icm. The temperature 
floor we apply (eq. 6) implies that the Field length is not de- 
fined for the cold phase in our simulations, and therefore that 
it is not relevant for convergence. Because the cold phase in our 
simulations does not cool, it can become pressure-supported 
at a finite size and resist further collapse. Convergence is some- 
what less restrictive in our simulations than in those studied 
by Sharma et al. (2010). 

We performed both 2d and 3d simulations and have 
confirmed that they give similar results. Many of the plots in 
this paper show the results of 2d simulations, since they are less 
expensive and permit a much larger parameter study. Because 
2d simulations contain fewer grid cells than 3d simulations, 
however, integrated quantities derived from 2d calculations 
are noisier. Thus, we chose to include only 3d simulations in 
Figs. 4 and 8. 

While our 2d and 3d simulations produce similar results, 
they are fundamentally different from one-dimensional simula- 
tions. Spatial variations between heating and cooling drive the 
local thermal instability in our model; hence, the development 
of multi-phase structure in our simulations is an inherently 
multi-dimensional effect. Additionally, the symmetry of a one- 
dimensional model prevents over-dense material from sinking 
and removes an important saturation channel from the ther- 
mal instability (§6). Much of the physics we describe in this 
paper is therefore absent in one-dimensional treatments of the 
icm such as those described in Ciotti & Ostriker (2001) and 
Guo & Oh (2008). 



5.4 Sensitivity to the Heating Function 

An important test of our model is the sensitivity of our con- 
clusions to the details of the (unknown) heating function. We 
study this dependence by adding random heating fluctuations 
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Figure 5. Convergence of the probability distribution function for density fluctuations P(5p/p). Colored lines show simulations at different 
resolutions, and the thick gray line shows the best-fit Gaussian distribution. This figure shows that the density inhomogeneities are reasonably 
converged in our simulations, apart from the obvious fact that one can resolve finer structure, and therefore higher Sp/p, at higher resolution. 
Note that in the limit of rapid cooling, the properties of the high density regions are determined in part by the temperature floor we apply, 
which determines the density of cold clumps that can be in pressure equilibrium with the surrounding hot plasma. 



of the form 

H^H(l + 5), (21) 

where S(x,t) is a Gaussian-random field with a white- noise 
spatial power spectrum and a temporal autocorrelation func- 
tion R$$(t) = e _r//tcorr . Thus, S introduces both spatial and 
temporal imbalances between heating and cooling, which per- 
sist for the coherence time ~ t CO rr- These fluctuations detune 
our feedback model while still preserving average thermal equi- 
librium, and are intended to mimic the temporal and spatial 
differences between heating and cooling which might arise in a 
more realistic feedback scenario. More importantly, including 
these fluctuations allows us to distinguish between results that 
are a consequence of the exact (and, in detail, unphysical) 
balance in equation 4 and results that are more robust and 
are primarily a consequence of global thermal stability. 

We carried out 2d simulations with £ cor r = {0.1, 1, 10}xtff 
and with the fluctuations normalized to root-mean-square 
(rms) amplitudes of 50%, 100% and 300%. (Note that we 
quote the RMS, or 'lcr' amplitude of the fluctuations; the peak 
values are considerably higher.) Fig. 6 shows images of the 
density fluctuations for these simulations and Fig. 7 shows the 
temperature distribution function for different values of t CO rr 
and the fluctuation amplitude. These figures demonstrate that 
our conclusions about the development of the thermal instabil- 
ity are essentially unaffected by order-unity fluctuations, over 
at least 10 cooling times. This important result implies that, 
as long as the plasma is in approximate global thermal equi- 
librium on reasonable time-scales ~ t Ti and length-scales ~ if, 



the development and saturation of local thermal instability 
will proceed approximately as shown in Figs. 1-4. We think 
that the existence of an approximate thermal equilibrium, 
rather than the specific details of our heating function (eq. 4), 
determines how the thermal instability develops and saturates. 
This conclusion is bolstered by Paper II, which finds very 
similar results using an entirely different heating function. 

Figs 6 and 7 show that extremely strong heating fluctu- 
ations with rms amplitudes of 300% spoil the thermal equi- 
librium of the plasma and induce a cooling catastrophe; even 
our extremely optimistic feedback model cannot withstand 
arbitrarily large heating perturbations. Though the feedback 
mechanism is not yet understood in clusters, this places a 
constraint on the heating: it should not differ persistently 
from the local cooling rate by more than a factor of several. 
Fig. 7 shows that this conclusion is essentially independent of 
the coherence time t CO rr of the heating. 



6 INTERPRETATION OF THE NON-LINEAR 
SATURATION 

In this section, we show that the linearized dynamical equa- 
tions provide valuable insight into the non-linear saturation of 
the thermal instability and its astrophysical implications. As 
in section 5, we focus on the development of multi-phase struc- 
ture (§6.2) and on the accreted mass flux (§6.3), which have 
been extensively studied observationally. Our basic procedure 
is to estimate a saturation amplitude for the linear instability. 
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Figure 6. Comparison of the plasma density in simulations with 
our fiducial heating function (equation 4) to simulations where we 
have added significant, random fluctuations to the heating function 
H (see eq. 21). These are white noise fluctuations with a temporal 
correlation t CO rr = fcrij simulations with longer correlation times 
^corr = 10£ti give similar results (see Fig. 7). From top to bottom, 
the panels show the density at t = 10 t TI in our fiducial simulations, 
simulations with 100% fluctuations in heating, and simulations 
with 300% fluctuations in heating. Fluctuations of 300% produce 
a cooling flow, but 100% fluctuations do not and instead produce 
results similar to our fiducial model. In all panels, color represents 
the log of the density, which ranges from 10~ 2 (blue) to 10 (red). 



Because we use linearized equations, the interpretation in this 
section only strictly holds in the weak cooling limit (£ T i ^> tff), 
so that the density perturbations remain relatively small. 

Fig. 8 illustrates the development and saturation of the 
thermal instability. The density inhomogeneity Sp/p (or any 
other quantity linear in the perturbation) initially grows ex- 
ponentially according to the dispersion relation (eq. 16), but 
eventually freezes out at a finite amplitude. This amplitude, 
along with the relations 12a-12c, then approximately deter- 
mines the state of the plasma at late times. In the following 
sections, we estimate this amplitude and show that we can re- 
produce elements of the non-linear saturation shown in Fig. 3. 
Fig. 8 shows that the difference between atmospheres which de- 
velop multi-phase gas and ones which do not is fundamentally 
a non-linear effect. The linear growth rate of the perturbations 
is largely independent of the time-scale ratio tn/ts, and it is 
the saturation which determines the degree of inhomogeneity 
at late times. 



6.1 Saturation Amplitudes 

In the limit that the plasma is buoyantly neutral (A = 0), we 
can estimate the saturation amplitude by inspecting the lin- 
earized, Lagrangian form of the momentum equation (eq. lb): 



The characteristic inflow (or outflow) time for a perturbed 
fluid element is t s i n k ~ H/Sv z . Initially, Sv z is small and this 
inflow time is long compared to the growth time of the thermal 
instability. As the perturbation grows, however, Sn/n increases 
and the fluid element accelerates according to equation 22. 
The inflow time £ S mk thus becomes shorter as the instability 
develops. We assume that the growth ceases when the inflow 
time is comparable to, or slightly shorter than, the growth time 
of the thermal instability. In this case, the thermal instability 
saturates when the velocity satisfies 

5v z ~!L. (23) 

6TI 

Thus, non-linear saturation occurs when a fluid element flows 
to smaller radii after one cooling time, as seems intuitively 
reasonable. 

The physical picture of a sinking fluid element does not 
apply in a stably stratified atmosphere, since the fluid element 
does not flow monotonically inwards, but instead oscillates 
with the gravity wave frequency A. The velocity associated 
with this oscillation dwarfs the mean, inward velocity. These 
waves are sourced by cooling, however, and we assume that 
they reach a steady-state in which the dissipation rate due 
to non-linear mode coupling equals the driving rate due to 
the thermal instability ~ t^ 1 . Thus the instability saturates 
when the dissipation time tdiss ~ H/Sv ~ t T i, where we have 
assumed strong turbulence and used the fact that the waves 
are driven on large scales, ~ H (as suggested by the bottom 
panels in Fig. 1). Though this saturation mechanism is very 
different from that described above for buoyantly neutral 
plasmas, it implies an equivalent saturation amplitude. We 
therefore assume that equation 23 describes the late-time 
evolution of the perturbations in all of our simulations. We 
show in the following sections how the behavior described 
in section 5 can be understood in terms of this saturation 
amplitude. 

6.2 Multi-Phase Structure 

Inserting our ansatz for the saturation amplitudes (equa- 
tion 23) into the momentum equation (12b) and using the 
dispersion relation (equation 16) to replace uj, we express the 
density inhomogeneity Sn/n at late times in terms of other 
properties of the plasma: 
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where tbuoy = A -1 and k]_ = (l — k 2 z ) is the squared horizontal 
component of the direction of the wave vector (typically ~ 1). 
Equation 24 has the asymptotic forms 
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^TI tbuoy 
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is 
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£ti ^> thu 



(25a) 
(25b) 



dv z 
~dt 
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Equation 25 shows that weakly stratified plasmas with t T i <C 
tbuoy develop smaller density inhomogeneities than plasmas 
with t T i ^> tbuoy (in the limit that t Tl > £ff). This difference 
arises because because plasmas with t T i ^> Woy can sustain 



(22) 



internal gravity waves, while atmospheres with t T 
cannot. 



< tbu 
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tn/ts = 10 tn/ts = 1 tn/ts = 1/10 




Figure 7. Gas mass as a function of temperature in simulations with different types of fluctuations about thermal equilibrium, measured at 
t = 10£ T i and z ~ H. These simulations are for isothermal initial conditions in which the initial temperature T = 1. The fluctuations are of 
the form H — > H(l + 5), where 5(x,t) has a white-noise spatial spectrum and a temporal coherence time t CO rr (see §5.4 for details). The gas 
properties are not sensitive to strong fluctuations in heating of up to 100% in amplitude. Stronger fluctuations of 300% generate significant 
cold material when t T i/% ^> 1; the cold material sinks to small radii (see Fig. 6), leading to a modest heating of the gas that remains at 
z ~ H. When £ T i/% <C 1, fluctuations of 300% break our ansatz of approximate thermal equilibrium and induce a cooling catastrophe. The 
correlation time £ C orr has only a modest influence on these results. 



Somewhat surprisingly, the right panel of Fig. 3 shows that 
the measured dependence o{5n/n on the time-scale ratio t T i/tff 
is in good agreement with equation 25b for both isothermal 
and isentropic initial conditions, even though tbuoy ^ oo in an 
isentropic atmosphere. This is because equation 25a applies 
only when the atmosphere is very nearly buoyantly neutral (at 
least when t T i is long compared to the dynamical time, as must 
be for equation 25 to be valid). While this is the case initially 
in our isentropic atmospheres, the entropy gradient evolves 
somewhat with time and equation 25a ceases to describe 
the plasma after only a few cooling times. The saturated 
values of tbuoy differ in our simulations with isentropic and 
isothermal initial conditions; equation 25b suggests that this 
may explain the systematic offset between these simulations 
shown in Fig. 3. For the longest cooling times, the evolution 
of the background profile is smallest; the slight steepening of 
Sp/p for our isentropic simulations in this limit may represent 
an intermediate case between equations 25a and 25b. 



6.3 Accreted Mass Flux 

We can also use our estimate of the non-linear saturation 
to understand the inward mass flux induced by the thermal 
instability. Defining the mass flux M — (5n5v z ) and using the 
estimates of Sv and Sn/n from the momentum equation and 



from equation 24, we find 
M —Mq¥ ( — V-4- 



x ( Re 



1 + 



tbu 



ik x \ t» / ik x\ 

e I Re ( e J 



(26) 



where Mcf = pH/t Tl is the mass flux expected in the absence 
of heating (recall that we are in the limit that t T i ^> Woy)- 
This yields 



M 
Mcf 



w 

















tff 

tbuoj 



(27) 



Equation 27 shows that the mass flux is dramatically sup- 
pressed when txi/tff ^> 1, qualitatively consistent with Fig. 4. 
Furthermore, this suppression is nearly independent of the ini- 
tial stratification of the plasma, even though stably stratified 
plasmas show much stronger density inhomogeneities (eq. 25). 
This is because the internal gravity waves that enhance 5n/n 
when £ T i/£buoy > 1 do not contribute to the net mass flux. 
Note that gravity waves dominate the instantaneous mass flux 
shown in Fig. 4 when t T i/tff ^> 1; the time-averaged accretion 
rate is smaller than the figure suggests. 
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Figure 8. Evolution of the density fluctuation Sp/p as a function 
of time in our simulations with isothermal initial conditions. The 
plotted quantity is an RMS average from z = 0.9H—1.1H. The 
density inhomogeneity grows from the initial perturbation until the 
characteristic infall time becomes comparable to the local cooling 
time. At this point, the density contrast saturates at approximately 
the value given by equation 25b. This figure shows the results from 
3d simulations, but the results from 2d simulations are similar. 



7 SIMULATIONS INCLUDING CONDUCTION 
AND MAGNETIC FIELDS 

The previous sections describe a simplified model of the ther- 
mal instability that neglects both conduction and the dynam- 
ical effect of the magnetic field. This model nicely isolates 
the physics of the thermal instability, but astrophysically it 
is too idealized. For example, conduction is critical for the 
thermal evolution of the plasma over a wide range of scales in 
the ICM, and Balbus (2000) and Quataert (2008) have shown 
that this completely changes the stability and dynamics of the 
plasma. In this section, we present results including magnetic 
fields and conduction and show that the conclusions from the 
previous sections largely apply in this more realistic case. 



7.1 Setup 

Our setup is very similar to that described in sections 2 and 3. 
We generalize equations lb and lc to include the effects of 
thermal conduction and the magnetic field: 



d_ 
dt 



. „ B z \ t B®B 
pv®v+ [P+—U+ — 

1 87T J 47T 



P9, 
(lb') 



pT 



ds 
dt 



(«-£)- V ■ Q cond , 

(lc') 



where B is the magnetic field and Q cond is the conductive 
heat flux. We evolve the magnetic field using the induction 



c)B 

dt 



= V x (v x B). 



(Id') 



We have ignored both (explicit) viscous and magnetic dissi- 
pation in equations lb'— Id'. These effects can influence MHD 
simulations in subtle and unexpected ways (see, e. g. Fromang 
& Papaloizou 2007; Davis et al. 2010), and so will need to be 
studied in detail in the future. 

The thermal conductivity of the plasma is strongly 
anisotropic in the icm and as a result the conductive heat flux 
is given by (Braginskii 1965) 



<3co„d = -nk BX eb(b-VT), 



(28) 



where b — B / B is a unit vector in the direction of the magnetic 
field and Xe is the thermal diffusivity of free electrons (with 
units of cm 2 /s). While the diffusivity Xe depends sensitively 
on temperature (Spitzer 1962), we take it to be constant in 
this exploratory analysis. This enables us to control the ratio 
of the conduction time to other timescales in the problem and 
thus to isolate the physics of cooling and conduction. Note 
that we still use the heating function denned by equation 4; 
any conductive heating or cooling of the plasma happens on 
top of the feedback heating. 

We initialize the plasma with a weak, horizontal magnetic 
field. (By 'weak,' we mean that magnetic tension is negligible 
in our simulations.) Because we impose reflecting boundary 
conditions at the upper and lower boundaries of the domain 
(§3), the magnetic field remains horizontal there and prohibits 
a conductive heat flux into the domain. 

As before, we solve equations la and lb'— Id' using 
Athena with the modifications described in section 3. We 
also implement equation 28 via operator splitting, using the 
anisotropic conduction algorithm described in Parrish & Stone 
(2005) and Sharma & Hammett (2007). In particular, we use 
the monotonized central difference limiter on transverse heat 
fluxes to ensure stability. This conduction algorithm is sub- 
cycled with respect to the main integrator with a time step 
At oc (Ax) 2 ; these simulations are therefore more computa- 
tionally expensive than adiabatic MHD calculations, especially 
at high resolution. 



7.2 Linear Properties 

We linearize equations la and lb'-ld' and perform a WKB 
analysis (see Quataert 2008 for more details). Assuming that 
magnetic tension is negligible, and proceeding as in section 4.1, 
the dispersion relation for the plasma is (cf. Balbus & Reynolds 
2010) 



p 3 - p 2 p F +pN 2 kl 



2 

^HBI 



0. 



(29) 



7 " 



In the above, p = —iuj is the growth rate of the perturbation, 
Pf = Pn — oj x is the growth rate of the thermal instability 
accounting for conduction (Field 1965), 

7 . . ; 2 (30) 

is inversely proportional to the conduction time across the 
wavelength of the perturbation, and 
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Figure 9. Comparison of the gas density in the non- linear state of simulations with different time-scale ratios £ti/% and £ x /£ff , including 
an initially horizontal magnetic field and anisotropic thermal conduction. We take the conduction time to be the time it takes heat to diffuse 
across one scale-height: t x = H 2 / x . Rapid conduction dramatically changes the morphology of the cold gas, smearing it out in the direction 
of the magnetic field. However, conduction does not appreciably change the mass of gas in the cold phase (Fig. 10). The ratio of timescales 
tn/ts still determines whether or not the plasma develops multi-phase structure, even in the limit of rapid conduction. In all frames, the 
color scale represents the log of the density; blue corresponds to a density of 10 -2 and red corresponds to a density of 10. The filaments 
are very straight in simulations where the Field length is longer than the domain size (right-most column) because conduction effectively 
eliminates all horizontal structure in the initial perturbations; the subsequent evolution is therefore nearly two-dimensional. We show both 
the top and bottom of the computational domain here to emphasize that the number of filaments/blobs produced by the thermal instability 
is somewhat stochastic. 



is the growth rate of either the magnet ot hernial instability 
(mti; Balbus 2000), or the heat-flux driven buoyancy instability 
(hbi; Quataert 2008). 

The MTI is unlikely to influence the development of multi- 
phase structure in galaxy clusters, since it operates outside 
the cool core, where the ratio of timescales £ T i/£ff is typically 
much greater than unity. While the hbi does operate efficiently 
in cool cores, it behaves like ordinary stable stratification in 
its saturated state (McCourt et al. 2011) and the growth time 
£ H bi is analogous to the timescale tbuoy used earlier. Thus, we 
do not expect the hbi or mti to change our results in any 
essential way (although this must be studied more carefully 
in future work). We anticipate that the same will be true for 
the overstabilities associated with the MTI and HBI (Balbus 
& Reynolds 2010). In this section, we use simulations with 
isothermal initial conditions (in which p H Bi — >• 0) so that these 
instabilities and overstabilities do not operate (at least in our 
initial conditions). This allows us to focus on the physics of 
thermal instability. 

The conduction frequency uj x is a function of scale, while 
the growth rate of the thermal instability p Ti is not; the modi- 



fied growth rate p F therefore must switch sign at the length- 
scale 



A F = \b • k\ x 



(2n) 



7 Pti 



1/2 



(32) 



known as the Field length (Field 1965). Intuitively, the Field 
length is the distance heat can diffuse in one cooling time; 
if the wavelength of a perturbation is larger than this dis- 
tance, conduction cannot stabilize it against cooling and the 
perturbation grows exponentially. 

Conduction suppresses the thermal instability on scales 
smaller than the Field length, but the Field length in a magne- 
tized medium depends on direction, as well as position. Even if 
the term in square brackets in equation 32 becomes arbitrarily 
large, the Field length will be small in directions orthogonal 
to the magnetic field. Because of this anisotropy, the thermal 
instability can still grow on scales much smaller than \/\etri- 
Sharma et al. (2010) have studied this growth in the absence of 
gravity; here we generalize their results to stratified plasmas. 

The growth rate for the thermal instability in our simula- 
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Figure 10. Mass fraction of cold material (with T < To/3) as a 
function of the timescale ratio tTi/% , for simulations with different 
conductivities. As in Fig. 3, this mass fraction is determined by 
averaging from z = 0.9-1.1 if and from t = 9-10 t TI . Lines show 
simulations with anisotropic thermal conduction and points indicate 
simulations with isotropic conduction. As suggested by Fig. 9, 
anisotropic thermal conduction does not strongly influence the 
multi-phase structure produced by the thermal instability: though 
the conductivities in these simulations differ by a factor of 100, 
the cold mass fractions agree to within about a factor of 2. In 
particular, the results appear to converge in the rapid-conduction 
limit. By contrast, we see no multi-phase structure in simulations 
with isotropic conduction if the Field length is comparable to, or 
greater than a scale-height (the arrows in this figure indicate upper 
limits on the mass fractions of cold gas). 



tions is 



\p F - N 2 kl/p F p F > \N\ 



-p F ±iNk± p F <^\N\ 



(33a) 
(33b) 



Equation 33 shows that the characteristic growth time of the 
thermal instability is Pp 1 , regardless of the entropy gradient. 
The growth rate p F reduces to p T i on large scales; thus, our 
results with and without conduction are very similar on scales 
larger than the Field length. Conduction prevents perturba- 
tions from growing below the Field length and therefore plays 
a similar role to the temperature floor in our non-conducting 
simulations. The primary difference between our conducting 
and non-conducting simulations is that the Field length is 
anisotropic in the conducting simulations, and the thermally 
unstable fluid elements collapse into long filaments, rather 
than the approximately spherical clumps shown in Figs. 1 
and 2. 



7.3 Numerical Results 

Fig. 9 shows 2d slices of 3d simulations with different values of 
the cooling constant Ao and the conductivity % e - We use only 
3d simulations in this section because, just as an over-dense 
fluid element cannot sink in one dimension, the dynamics of 
a sinking magnetized filament changes in going from two to 
three dimensions. These simulations all use our isothermal 
initial condition and initially have weak, horizontal magnetic 
field lines in the plane of the figure. Rapid conduction smears 
out the cold clumps into filaments of length ~ Af, but does not 
otherwise alter the growth of the thermal instability. Specifi- 
cally, Fig. 9 demonstrates that, even in the limit of very rapid 
conduction, the ratio of time-scales t TI / tff determines whether 
the plasma develops multi-phase structure. This result depends 
critically on the anisotropic nature of thermal conduction. In 
the rightmost panels of Fig. 9, the Field length is larger than 
the entire simulation domain; if conduction were isotropic, the 
entire atmosphere would become nearly isothermal and the 
thermal instability would be suppressed. The insulating ef- 
fect of the magnetic field permits large temperature gradients 
orthogonal to the magnetic field and thus the formation of 
multi-phase structure (Sharma et al. 2010). 

Fig. 10 quantifies the effect of conduction on the thermal 
instability: we show the mass fraction of cold gas (as in the 
left panel of Fig. 3) for 3d simulations with different ther- 
mal conductivities. In simulations with anisotropic thermal 
conduction, this mass fraction is almost independent of the 
conductivity, and it appears to converge in the limit that the 
conductivity becomes large. This behavior is consistent with 
Fig. 9. Together, these results imply that anisotropic conduc- 
tion alters the morphology of the gas in the cold phase, but 
not the presence, absence, or amount of multi-phase structure. 

We have also run a number of simulations with isotropic 
thermal conduction. These simulations use the same setup as 
before, but with the conductive heat flux Q cond — — nfeXe VT, 
where (as before) % e is a constant, free parameter. In order to 
prevent conduction from changing the total energy content of 
the plasma, we set %e = at the upper and lower boundaries of 
the computational domain so that there is no conductive heat 
flux into the domain. Fig. 10 shows that, while anisotropic 
conduction does not strongly influence the amount of cold 
gas produced by the thermal instability, isotropic conduction 
can quench it entirely: we see no multi-phase structure in 
our simulations with isotropic conduction whenever the Field 
length is comparable to, or larger than, the pressure scale- 
height. These conclusions also apply to other properties of the 
plasma quantified in section 5, e.g. the accreted mass flux: 
anisotropic thermal conduction has little effect on this quantity, 
while isotropic thermal conduction can strongly suppress it. 

Voit et al. (2008) suggested that thermal instability pro- 
duces multi-phase gas in clusters when the Field length is 
comparable to, or smaller than, the size of the cool core, 
but that conduction suppresses the formation of multi-phase 
structure for larger Field lengths. Coincident ally, in typical 
cool-core clusters, this criterion is quantitatively similar to 
our criterion on the ratio t T i/tff. 3 However, because the icm 
is magnetized, thermal conduction is extremely anisotropic; 



3 This comparison makes use of the result from Paper II that 
the threshold for multi-phase gas in spherical systems is closer to 
tn/tff ~ 10. 
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Table 2. Clusters used in Fig. 11. 



Extended H Q 


No Extended H a 


Abell 133 


Abell 85 


Abell 478 


Abell 644 


Abell 496 


Abell 744 


Abell 780 


Abell 1650 


Abell 1795 


Abell 2029 


Abell 1991 


Abell 2142 


Abell 2597 


Abell 4059 


Sersic 159-03 




Centaurus 





We use the surveys of McDonald et al. (2010, 2011a,b) to determine 
whether a cluster shows multi-phase gas, and we use the data in 
the ACCEPT catalog to estimate t T i/% for the hot ICM. Our label 
'extended H a ' signifies that the H Q emission can be resolved and is 
known to exist outside the BCG; these are the Type I systems from 
McDonald et al. 

the results of this section demonstrate that even very rapid 
thermal conduction cannot suppress local thermal instabil- 
ity. Thermal conduction only stabilizes modes parallel to the 
magnetic field, and multi-phase structure continues to develop 
via perturbations that are roughly orthogonal to the local 
magnetic field. 



8 DISCUSSION 

Observational limits from x-ray spectroscopy (Peterson & 
Fabian 2006) and from the shape of the galaxy luminosity 
function (Benson et al. 2003) indicate that the diffuse plasma 
in galaxy groups and clusters does not cool as quickly as it ra- 
diates. These observations imply that some heating process off- 
sets radiative cooling and that the gas remains in approximate 
thermal equilibrium, at least when averaged over length-scales 
comparable to the scale-height or time-scales comparable to 
the cooling time. The nature of this "feedback" heating is not 
yet fully understood, although it appears to involve heating by 
a central AGN (e. g. Birzan et al. 2004; McNamara & Nulsen 
2007). More fully understanding the mechanism(s) that regu- 
late the heating to so closely match cooling remains a major 
challenge in theories of galaxy formation. 

Though heating strongly suppresses cooling in galaxy 
groups and clusters, star formation and multi-phase gas pro- 
vide clear evidence for cold gas in many cluster cores. Obser- 
vational indicators of this cold gas strongly correlate with the 
cooling time of the ambient hot ICM (e. g. Voit et al. 2008; 
Rafferty et al. 2008; Cavagnolo et al. 2008, 2009), motivating 
a model in which thermal instability in the hot ICM produces 
much of the cold gas in cluster cores (as has been suggested 
many times in the past, e.g. Fabian & Nulsen 1977; Cowie 
et al. 1980; Nulsen 1986; Loewenstein et al. 1991). Theoreti- 
cally studying local thermal instability in the icm has proven 
difficult, however, because of the cooling-flow problem: studies 
that include both cooling and gravity typically find that the 
plasma is globally thermally unstable, and that the entire 
cluster core collapses monolithically. This difficulty has led 
several authors to conclude that the thermal instability does 
not produce multi-phase structure in stably-stratified systems 
at all. 

We avoid the cooling-flow problem in this paper by adopt- 



ing a phenomenological heating model that enforces thermal 
equilibrium when averaged over large scales (§2.1). Our heat- 
ing model is approximate, over-simplified, and wrong in detail. 
However, our results are insensitive to large temporal and spa- 
tial fluctuations about the average heating (§5.4). Moreover, 
in Paper II we obtain similar results using a more physical 
feedback heating prescription. We therefore believe that our 
conclusions about the saturation of the local thermal instabil- 
ity are reasonably robust. 

In the current paradigm in which clusters are approxi- 
mately in global thermal equilibrium, heating of the ICM is 
very likely to depend explicitly on position in the cluster. In 
this case, the thermal stability of the plasma is independent of 
its convective stability (§4.2). Fundamentally, buoyancy drives 
convection, while heating and cooling drive thermal instabil- 
ity. These two processes are formally related only under the 
restrictive assumption that heating is a state function of the 
plasma; more generally, they are unrelated and the thermal 
stability of a plasma is independent of its convective stability 
(see Balbus & Soker 1989). 

Nonetheless, it remains true that the competition between 
buoyancy and thermal instability determines the net effect of 
cooling on a stratified plasma. We parametrize the relative 
importance of these effects using the dimensionless ratio t T i/£ff 
(the ratio of the thermal instability growth time txi to the 
local dynamical, or free-fall, time tff). When this ratio is 
small, thermal instability dominates and the plasma develops 
significant multi-phase structure; when the ratio £ T i/£ff is large, 
buoyancy dominates and the plasma remains in a single, hot 
phase (§5). This dependence of the saturation on £ T i/£ff is 
true for both stably stratified and neutrally stratified plasmas 
(Fig. 3), even though the effect of buoyancy is very different 
in these two cases. 

More quantitatively, we find (Fig. 3) that the saturated 
density inhomogeneities produced by the thermal instability 
approximately obey the relation 

this scaling can be understood analytically by assuming a 
saturation amplitude for the thermal instability in which the 
characteristic fluid velocities approach v sa t ~ H/t coo \ (§6). 
Thus, by assuming that some heating mechanism prevents 
cooling catastrophes in clusters, we find that the ICM breaks up 
into multiple phases via local thermal instability (with Sp/p > 
1) only if the dimensionless ratio of timescales £ T i/£ff < 1; 
specifically, there is almost no cold gas at large radii when 
£ti > tff (Fig. 3). This finding is one of the primary results of 
our analysis. We note that the linear growth of the thermal 
instability is largely independent of the ratio t T i/tff. Thus the 
difference between atmospheres which develop multi-phase 
gas and those which do not is fundamentally due to how the 
non-linear saturation of the thermal instability depends on 
the atmosphere's properties. 

The calculations in Paper II show that the criterion for 
multi-phase structure is actually somewhat less stringent in 
spherical systems, £ T i/£ff < 10. This difference stems from 
the fact that fluid elements are compressed as they move 
inwards in a spherical system; this compression enhances 
the density perturbations and accelerates the growth of the 
thermal instability. 

Our criterion for multi-phase structure is not sensitive 
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Figure 11. {Left:) The timescale ratio t C oolAff as a function of radius for clusters in both the ACCEPT catalog and the McDonald et al. 
(2010) survey. Solid blue lines show clusters with filaments and dashed red lines show clusters that lack detected extended H a emission. 
Clusters with filaments have systematically lower values of t C oolAff- Furthermore, this ratio is smallest between ~10— 50 kpc, where most 
filaments are found. (Right:) The same clusters in the t C ool _ tcoolAff plane. The coloring is the same as in the left panel. The ratio £ C oolAff 
appears to be a slightly better predictor of multi-phase structure than t coo \ alone. Table 2 lists the clusters plotted in this figure. 



to large variations about our idealized heating prescription 
(§5.4 and Fig. 6). Furthermore, the multi-phase structure 
that develops via thermal instability is largely independent 
of the magnitude of the thermal conduction, even on scales 
much smaller than the Field length (provided conduction is 
anisotropic, as is the case in galaxy groups and clusters; §7). 
Anisotropic thermal conduction changes the morphology of the 
cold gas produced via thermal instability (blobs — >> filaments), 
but not the presence or amount of cold gas. There is thus 
a very strong, qualitative difference between isotropic and 
anisotropic thermal conduction (Fig. 10), which cannot be 
captured by simply multiplying the heat flux by a suppression 
factor (as is often done, e. g. Zakamska & Narayan 2003; Voit 
et al. 2008; Guo & Oh 2008). 

Our heating prescription imposes global thermal equilib- 
rium and reduces the accreted mass flux in our model halos 
relative to cooling-flow values. This reduction is not inevitable, 
however, because thermal instability can produce cooling-flow- 
like inflow rates when t TI < tff. In a globally stable system, the 
thermal instability thus plays an important role in regulating 
gas inflow rates. We study the connections among thermal 
instability, mass inflow and feedback more fully in Paper II. 

We argue that a locally stable heating mechanism (such 
as the one proposed in Kunz et al. 2011) is not required to ex- 
plain the reduced star formation and cooling rates in clusters. 
Instead, global stability arising from approximate thermal 
equilibrium, together with the physics of local thermal insta- 
bility in stratified plasmas, is sufficient to reproduce the low 
net cooling rates in clusters. Moreover, the correlation of H a 
filaments and star formation in clusters with the cooling time 
in the hot icm strongly suggests that the plasma is in fact 
locally thermally unstable. We have shown that the suppres- 
sion of accretion rates is not sensitive to thermal conduction 



or to significant variations about our specific feedback pre- 
scription (§5.4 and Fig. 6), and we explore this further in 
Paper II. 

We now compare our model predictions with observational 
results. The thermal instability time t T i = Pti 1 (eq. 19) depends 
on the unknown way in which feedback energy is thermalized 
(eq. 4) and cannot be directly inferred from observations. We 
therefore use the cooling time when comparing our results 
with observations. These two timescales differ by an unknown 
factor of order unity. The ratio of timescales t C ooi/£ff can be 
reexpressed in more familiar terms using 



i cool (if/10 keV cur 



2\3/2 



T 7 1/2 A- 23 Myr)' 



(35) 



where K = /cbT/tt 2 / 3 is a measure of the plasma entropy, T7 
is the temperature in units of 10 7 K and A_23 is the cooling 
function (eq. 5) in units of 10 -23 erg cm 3 s _1 . Thus, the 
plasma in clusters and galactic halos should show extended 
multi-phase structure wherever 



K < (20 keV cm 2 



-a/2 



tff 



30 Myr 



2/3 



(36) 



where we have used the threshold from Paper II for multi-phase 
structure in spherical systems: £ T i/£ff ~ 10. This criterion is 
consistent with observations that clusters with central en- 
tropies below 30 keV cm 2 preferentially show signs of cold gas 
such as star formation and H a emission (Voit et al. 2008; Cav- 
agnolo et al. 2008, 2009). Note, however, the relatively strong 
dependence of this criterion on tff and on A; this is because 
the entropy K is not the fundamental parameter governing 
the saturation of the thermal instability in a stratified system. 
The Hq, survey conducted by McDonald et al. (2010, 
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2011a) permits another test of our criterion. McDonald etal. 
provide lists of groups and clusters with and without extended 
Hq, emission. We test our criterion by estimating the time- 
scale ratio tcooi/tff for these systems using data from the 
accept catalog (Cavagnolo et al. 2009). We fit the entropy 
profiles of clusters in the ACCEPT catalog using K(r) = Kq + 
i^i(r/100kpc) a and we fit the pressure profiles P(r) using 
the form provided in Arnaud et al. (2010). Of the groups and 
clusters in both the H a surveys and in accept, sixteen give 
reasonable fits (listed in Table 2). 4 From our fits to K(r) and 
P(r), we calculate n(r) and T(r) and estimate t coo \(r) using 
the fit to the cooling function provided by Tozzi & Norman 
(2001) with 1/3 solar metallicity. We estimate g(r) from the 
pressure and density profiles by assuming spherical symmetry 
and hydrostatic equilibrium; from this we calculate tff(r). 

The left panel of Fig. 11 shows t coo i/tff as a function of 
radius for these sixteen groups and clusters. As predicted by 
our analysis, the clusters with short cooling times £ CO oi/£ff ^ 10 
show extended filaments, while the clusters with long cooling 
times £ C ooiAff ^ 10 do not. Additionally, most of the filaments 
are found at radii around 10-50 kpc, where the ratio £ T i/£ff is 
the smallest. Although this evidence is not conclusive, these 
data support our hypothesis that the filaments condense from 
the icm due to the local thermal instability. 

We emphasize that both cooling and gravity influence the 
development of the thermal instability in the icm. A short 
cooling time (or low K) is not sufficient for the formation of 
filaments; rather, the ratio £ T i/£ff is the relevant parameter. 
The right panel of Fig. 11 shows clusters in the (t CO oi)-(t C ooi/tff) 
plane. More data are needed to conclusively test our model, 
but these results are consistent with our interpretation that the 
ratio of timescales t CO oi/£ff is a better predictor of multi-phase 
gas in hot halos than t coo i alone. 

These simple comparisons support a model in which local 
thermal instability produces at least some of the H a filaments 
seen in clusters. This does not, however, imply that thermal 
instability alone can explain all of the observed properties of 
multi-phase gas in clusters and/or galaxies. On the contrary, 
processes such as conductive condensation of hot gas to cool 
gas ("non-radiative cooling;" Fabian et al. 2002; Soker et al. 
2004) and the inflow of cold gas through the virial radius (Keres 
& Hernquist 2009) may also be important (in higher and lower 
mass halos, respectively). Furthermore, other processes in the 
icm such as merger shocks, galaxy wakes and buoyant radio 
bubbles may influence the evolution of the filaments. 

Accretion of the cold gas formed via thermal instability 
likely plays an important role in the evolution of brightest 
cluster galaxies and their central AGN (Pizzolato & Soker 2005, 
2010). In addition, the high- velocity clouds surrounding the 
Milky Way may also be manifestations of the thermal instabil- 
ity (Mailer & Bullock 2004; Sommer-Larsen 2006; Kaufmann 
et al. 2006; Peek et al. 2008); this process could provide an 
important source of unenriched gas to maintain metallicity 
gradients (Jones et al. 2010) and continued star formation 
(Bauermeister et al. 2010) in the Milky Way and other galaxies. 

4 Unfortunately, we were unable to fit several clusters with well- 
known filament systems, including Perseus, Abell 2052 and M87. 
The pressure gradients in Perseus and Abell 2052 are positive at 
some radii and cannot be fit by the universal profile. Similarly, 
the pressure profile for M87 deviates from the broken power-law 
universal profile. 
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